# Replicates all analyses in Face-to-face vs. Virtual Assistance to Formalize: Evidence From a Field Experiment in Brazil\
# By Zucco, Lenz, Goldszmidt, and Valivia paper submitted to Economic Letters
# Generates all tables and figures in main paper and in supplemental information

library(haven)
library(tidyverse)
library(plyr)
library(data.table)
library(DeclareDesign)
library(estimatr)
library(randomizr)
library(xtable)
setwd("~/Dropbox/Data/Metaketa/")

load(file="ProcessedData/Data_Letters.RData")
ed <- subset(dbrazil,Z==0|Z==1|Z==2)
ed$Z <- factor(ed$Z)

# Differentiate between compliance with T1 and with T2
ed$D[which(ed$Z==2&ed$D2==1)] <- 2
ed$D <- factor(ed$D)
print(table(D=ed$D,Z=ed$Z))

# Write in probabilities of assignment (stable across all blocks)
ed$Z_prob <- as.numeric(car::recode(ed$Z,"0=.375;1=.375;5=.250"))

## Analysis ##

## Declare functions used in the analysis ##
my.summary <- function(index, the.data=ed){
  out <- as.data.frame(matrix(NA, 1, 6))
  tmp <- summary(the.data[,outcomes[index]],na.rm=T)
  out[1,1] <- tmp["Mean"]
  out[1,2] <- sd(the.data[,outcomes[index]],na.rm=T)
  out[1,3] <- tmp["Min."]
  out[1,4] <- tmp["Max."]
  out[1,5] <- tmp["NA's"]
  out[1,6] <- length(table(the.data[,outcomes[index]]))
  names(out) <- c("Mean","SD","Min","Max","Missing","Type") 
  return(out)
}

## These are set up simply to spead up estimation of each model for all outcome variables 
est_dm <- function(index, treat = "Z", the.data=ed){
  require("ANOVAreplication")
  results_t <- as.data.frame(matrix(NA, 1, 16))
  the.formula<-formula(paste(outcomes[index],"~Z"))
  dm1 <- difference_in_means(the.formula, data =  the.data, condition1=0,condition2=1)
  dm2 <- difference_in_means(the.formula, data =  the.data, condition1=0,condition2=2)
  dm21 <- difference_in_means(the.formula, data =  the.data, condition1=1,condition2=2)
  results_t[1:7] <- summary(dm1)$coef
  results_t[8:14] <- summary(dm2)$coef
  results_t[15] <- summary(dm21)$coef[,4]
  results_t[16] <- dm1$N
  names(results_t) <- c(paste(colnames(summary(dm1)$coef),1,sep="."),
                        paste(colnames(summary(dm2)$coef),2,sep="."),"p.T2vT1","N")
  cat(outcomes[index],"\n")
  return(results_t)
}

#Similar, but this variant is for attrition balance
est_attr <- function(index, treat = "attrition",the.data=ed,w=NULL,std=F){
  require("ANOVAreplication")
  results_t <- as.data.frame(matrix(NA, 1, 8))
  the.formula<-formula(paste(outcomes[index],"~",treat))
  dm <- difference_in_means(the.formula, data = the.data,alpha = 0.1,weights=w)
  results_t[1:7] <- summary(dm)$coef
  results_t[8] <- dm$N
  names(results_t) <- c(colnames(summary(dm)$coef),"N")
  if(std==T){#compute standardized mean differences
    tmp <- data.frame(y = the.data[,outcomes[index]], g = the.data[,treat])
    pool_sd <- pooled.sd(na.omit(tmp))
    results_t <- cbind(results_t,'Std. Estimate'=round(results_t[1,1]/pool_sd,5))
  }
  return(results_t)
}

# Basic ITT with Z0 compared to Z1 and Z2
est_itt <- function(index, treat = "Z", ctrls=NULL, fe=F, baseline=F,the.data=ed){
  results_t <- as.data.frame(matrix(NA, 1, 15))
  if(is.null(ctrls)){ 
    the.formula <- formula(paste(outcomes[index]," ~ Z"))
  }else{
    if(baseline==F){ the.formula <- formula(paste(outcomes[index],"~ Z +",paste(ctrls,collapse="+")))}
    if(baseline==T){ the.formula <- formula(paste(outcomes[index],"~ Z +",paste(ctrls,collapse="+"),"+",baselines[index])) }
  }
  if(fe){  
    reg <- lm_robust(the.formula, se_type="stata", fixed_effects= ~ block, data = the.data)
  }else{reg <- lm_robust(the.formula, data = the.data, se_type="stata")
  }
  results_t[1:7] <- round(summary(reg)$coef["Z1",],5)
  results_t[8:14] <- round(summary(reg)$coef["Z2",],5)
  results_t[15] <- reg$N
  names(results_t) <- c(paste(colnames(summary(reg)$coef),1,sep="."),
                        paste(colnames(summary(reg)$coef),2,sep="."),"N")
  cat(outcomes[index],"\n")
  return(results_t)
}

# Heterogenous effects ITT with Z2 unfolded in to Z2l and Z2h
est_itth <- function(index, treat = "Zi", ctrls=NULL, fe=F, baseline=F,the.data=ed){
  results_t <- as.data.frame(matrix(NA, 1, 22))
  if(is.null(ctrls)){ 
    the.formula <- formula(paste(outcomes[index]," ~ Zi"))
  }else{
    if(baseline==F){ the.formula <- formula(paste(outcomes[index],"~ Zi +",paste(ctrls,collapse="+")))}
    if(baseline==T){ the.formula <- formula(paste(outcomes[index],"~ Zi +",paste(ctrls,collapse="+"),"+",baselines[index])) }
  }
  if(fe){  
    reg <- lm_robust(the.formula, se_type="stata", fixed_effects= ~ block, data = the.data)
  }else{reg <- lm_robust(the.formula, data = the.data, se_type="stata")
  }
  results_t[1:7] <- round(summary(reg)$coef["Zi1",],5)
  results_t[8:14] <- round(summary(reg)$coef["Zi2",],5)
  results_t[15:21] <- round(summary(reg)$coef["Zi3",],5)
  results_t[22] <- reg$N
  names(results_t) <- c(paste(colnames(summary(reg)$coef),1,sep="."),
                        paste(colnames(summary(reg)$coef),2,sep="."),
                        paste(colnames(summary(reg)$coef),3,sep="."),"N")
  cat(outcomes[index],"\n")
  return(results_t)
}

est_cace <- function(index, treat = "Z", ctrls=NULL, fe=F, baseline=F, comply="D",the.data=ed){
  the.data$D <- the.data[,comply] 
  if(index==1){cat("Separate IV estimation for Z1 and Z2 relative to Z0!\nUsing"
                   ,comply,"as indicator of compliance\n")}
  results_t <- as.data.frame(matrix(NA, 1, 15))
  if(is.null(ctrls)){ 
    the.formula <- formula(paste(outcomes[index]," ~ D  | Z"))
  }else{
    if(baseline==F){ the.formula <- formula(paste(outcomes[index],"~ D +",paste(ctrls,collapse="+"),"| Z +",paste(ctrls,collapse="+")))}
    if(baseline==T){ the.formula <- formula(paste(outcomes[index],"~ D +",paste(ctrls,collapse="+"),"+",baselines[index],"| Z +",paste(ctrls,collapse="+"),"+",baselines[index])) 
    }
  }
  if(fe){  
    regZ1 <- iv_robust(the.formula, se_type="stata", fixed_effects= block, data = subset(the.data,Z==0|Z==1))
    regZ2 <- iv_robust(the.formula, se_type="stata", fixed_effects= block, data = subset(the.data,Z==0|Z==2))
  }else{regZ1 <- iv_robust(the.formula, data = subset(the.data,Z==0|Z==1), se_type="stata")
  regZ2 <- iv_robust(the.formula, data = subset(the.data,Z==0|Z==2), se_type="stata")
  }
  results_t[1:7] <- round(summary(regZ1)$coef["D1",],5)
  results_t[8:14] <- round(summary(regZ2)$coef["D2",],5)
  results_t[15] <- regZ2$N
  names(results_t) <- c(paste(colnames(summary(regZ1)$coef),"1",sep="."),
                        paste(colnames(summary(regZ2)$coef),"2",sep="."),"N.2")
  cat(outcomes[index],"\n")
  return(results_t)
}

# Pre-Analysis Recodings ######################################
# Replace -999,-998, and -997 values with NA's
ed[ed == -997|ed == -998|ed==-999] <- NA

# Rescale income to facilitate presentation
ed$income_survey_b <- log(1+ed$income_survey_b)
ed$income_survey_e <- log(1+ed$income_survey_e)

# Rescale know_meifeatures
ed$know_meifeaturesi_e <- ed$know_meifeatures_e
ed$know_meifeatures_e <- as.numeric(ed$know_meifeatures_e > mean(ed$know_meifeatures_e[ed$Z==0],na.rm=T))

# Rocode tax compliance
ed$taxcompliance_admin_e[which(ed$formalize_admin_e==0)] <- 0

# Declare control variables #
the.ctrls <- c("incomebracket_survey_b"#income brackets
               ,"female_survey_b"#sex
               ,"white_b"#race
               ,"bussize_b"#size of business
               ,"bussector_b"#sector of activity
               ,"bustype_b"#type of business (fixed store or street vendor)
) 

# Examine missingness in control variables (on non-attrited cases)
as.matrix(apply(is.na(subset(ed
                             ,select=the.ctrls,is.na(intenttoformalize_admin_e)==F|
                               is.na(formalize_admin_e)==F|
                               is.na(know_meigeneral_e)==F)),2,sum))

# Create additional "missing" category for NA observatins in control variables with > 10 NAs
# Make sure controls with more than two categories are treated as categorical variables
ed$incomebracket_survey_b[which(is.na(ed$incomebracket_survey_b))] <- -999
ed$incomebracket_survey_b <- factor(ed$incomebracket_survey_b)
ed$white_b[which(is.na(ed$white_b))] <- -999
ed$white_b <- factor(ed$white_b)

# Attrition #################################################

# Compare attrition indicators
table(Survey=ed$attrition,Admin=ed$attrition_admin)

# Compliance among non-attriters
table(D=ed$D,SurveyAttr=ed$attrition,Z=ed$Z)

# Compliance among non-attriters
table(D=ed$D,AdminAttr=ed$attrition_admin,Z=ed$Z)

# Examine attrition (endline survey)
reg_at1 <- lm_robust(attrition~Z,se_type="stata",data=ed)  
print(reg_at1)

# Examine attrition (administrative data)
reg_at2 <- lm_robust(attrition_admin~Z,se_type="stata",data=ed)  
print(reg_at2) 

# Examine attrition (administrative data+direct observation)
reg_at3 <- lm_robust(is.na(intenttoformalize_admin_e)~Z,se_type="stata",data=ed)  
print(reg_at3) 

#Plot differential attrition
png(file = "Figures/figure-letters-differentialattrition.png",
    width = 7, height = 8, units= "in" ,res = 300)
par(mar=c(3,4,2,.5))
x1 <- summary(reg_at1)$coef
x2 <- summary(reg_at2)$coef
to.plot <- cbind(c(x1[1,1]+x1[2,1],x1[1,1]+x1[3,1])
                 ,c(x2[1,1]+x2[2,1],x2[1,1]+x2[3,1]))
tmp <- barplot(to.plot,ylim=c(0,.4),beside=T
               ,ylab="Share of Sample"
               ,legend.text=c("In Person","Inst. Msg.")
               ,args.legend = list(x="topright",bty="n")) 
segments(x0=tmp[1,1]-0.6, x1=tmp[2,1]+0.6, y0=x1[1,1],y1=x1[1,1],lty=2)
segments(x0=tmp[1,2]-0.6, x1=tmp[2,2]+0.6, y0=x2[1,1],y1=x2[1,1],lty=2)
axis(side=1,at=apply(tmp,2,mean),labels=c("Endline","Administrative"),tick=F)
dev.off()

# Examine balance between attrited and non-attried on baseline variables
aoutcomes <- c("access_survey_b",
               "access_publicse_survey_b",
               "rses_b",
               "life_b",
               "intpertrust_b",
               "entrepreneurship_b",
               "managerialstyle_b",
               "female_survey_b",
               "age_survey_b",
               "education_survey_b",
               "income_survey_b")
ann <- length(aoutcomes)    

#F-test predicting attrition
the.formula <- formula(paste("attrition ~",paste(aoutcomes ,collapse="+")))
summary(lm(the.formula,data=ed))
the.formula <- formula(paste("attrition_admin ~",paste(aoutcomes ,collapse="+")))
summary(lm(the.formula,data=ed))

#Compute difference in means for each baseline variable on survey attrition
outcomes <- aoutcomes 
attr_1 <- bind_rows(lapply(1:ann, FUN = est_attr, treat="attrition",std=T))
rownames(attr_1)<-outcomes
#print(round(attr_2[,-c(5:7)],4))

#Compute difference in means for each baseline variable on admin attrition
attr_2 <- bind_rows(lapply(1:ann, FUN = est_attr, treat="attrition_admin",std=T))
rownames(attr_2)<-outcomes
#print(round(attr_2[,-c(5:7)],4))

# Plot the same information
png(file = "Figures/figure-letters-attritionbalance.png",
    width = 7, height = 8, units= "in" ,res = 300)
par(mar=c(4,6,3,.5))
bal1 <-attr_1[,"Std. Estimate"]
bal2 <-attr_2[,"Std. Estimate"]
plot(bal1,-1*(1:nrow(attr_1)+0.05),bty="n"
     ,xlab="Standardized Mean Differences",ylab="",xlim=c(-.6,.6)
     ,yaxt="n",  pch=ifelse(attr_1[,4]>0.1,22,24),xpd=NA,main="Balance on Attrition")
points(bal2,-1*(1:nrow(attr_1)-0.05), pch=ifelse(attr_1[,4]>0.1,22,24),bg=gray(0.5))
abline(v=c(-.25,0,.25),lty=c(2,1,2),col=c(gray(0.5),1,gray(0.5)))
axis(side=2,at=-1*1:nrow(attr_1),labels=rownames(attr_1),tick=F,las=2,line=-2,xpd=NA,cex.axis=0.7)
legend(x="right",legend=c("Endline","Admin"),pt.bg=c(0,gray(0.5)),pch=c(22,22),bty="n",cex=0.7)
dev.off()

### Examine ineligibility
print(prop.table(table(Ineligible=ed$ineligible)))
print(table(Ineligible=ed$ineligible,Attrition=ed$attrition))
reg_i <- lm_robust(ineligible~Z,se_type="stata",data=ed) 
print(reg_i)

the.formula <- formula(paste("ineligible ~",paste(aoutcomes ,collapse="+")))
summary(lm(the.formula,data=ed))

# Main outcomes used in Econ Letters Endorsement paper
paper.outcomes <- c("know_meigeneral_e",
                    "know_meifeatures_e",
                    "intenttoformalize_admin_e",
                    "formalize_admin_e",
                    "formalizeall_admin_e",
                    "taxcompliance_admin_e" 
					)

all.pbaselines <- c(1,1,"intenttoformalize_survey_b",1,1,1)
pn <- length(paper.outcomes)
outcome.labels <-
  c("Knoweldge (General)"
    ,"Knowledge (Features)"
    ,"Intent to Formalize" 
    ,"Formalize (as MEI)"
    ,"Formalize (any)"
    ,"Formalize & Comply")

### Create table of descriptives statisticas for paper
outcomes <- c(paper.outcomes)
the.summary <- bind_rows(lapply(1:length(outcomes), FUN = my.summary, the.data=ed))
rownames(the.summary) <- outcome.labels
the.summaryx <- xtable(the.summary, digits=c(0,2,2,0,0,0,0))
caption(the.summaryx) <- "Descriptive Statistics" 
label(the.summaryx) <- "tab:descriptive" 
print(the.summaryx,caption.placement = "top")

# Estimation using only variables used in paper, but otherwise following PAP #
# Raw Diferences in Outcome # 
outcomes <- paper.outcomes
res_1 <- bind_rows(lapply(1:length(outcomes), FUN = est_dm, the.data=ed))
rownames(res_1) <-  outcomes
print(round(res_1[,c(1,4,8,11,15,16)],4))

# ITT #
# All outcomes #
outcomes <- paper.outcomes
nn <- length(outcomes)
baselines <- all.pbaselines
res_3a <-  bind_rows(lapply(1:nn, FUN = est_itt,fe=T)) 
res_3b <-  bind_rows(lapply(1:nn, FUN = est_itt,fe=F,ctrls=the.ctrls))  
res_3c <-  bind_rows(lapply(1:nn, FUN = est_itt,fe=T,ctrls=the.ctrls,baseline=T))  
res_3d <-  bind_rows(lapply(1:nn, FUN = est_itt,fe=F
                            ,ctrls=c(the.ctrls,"ineligible"),the.data=ed)) 
res_3e <-  bind_rows(lapply(1:nn, FUN = est_itt,fe=T
                            ,ctrls=c(the.ctrls,"ineligible"),the.data=ed)) 
rownames(res_3a) <- rownames(res_3b) <- rownames(res_3c) <-  
  rownames(res_3d) <- rownames(res_3e) <- outcomes 
print(round(res_3d[,c(1:4,8:11)],4))

#  CACE (omitted from PAP)
outcomes <- paper.outcomes
nn <- length(outcomes)
baselines <- all.pbaselines
res_5a <-  bind_rows(lapply(1:nn, FUN = est_cace,fe=T,the.data=ed))  
res_5b <-  bind_rows(lapply(1:nn, FUN = est_cace,fe=F,ctrls=the.ctrls,the.data=ed))
res_5c <-  bind_rows(lapply(1:nn, FUN = est_cace,fe=T,ctrls=the.ctrls,baseline=T,the.data=ed)) 
res_5d <-  bind_rows(lapply(1:nn, FUN = est_cace,fe=F
                            ,ctrls=c(the.ctrls,"ineligible"),baseline=T,the.data=ed)) 
res_5e <-  bind_rows(lapply(1:nn, FUN = est_cace,fe=T
                            ,ctrls=c(the.ctrls,"ineligible"),baseline=T,the.data=ed)) 
rownames(res_5a) <- rownames(res_5b) <- rownames(res_5c) <- rownames(res_5d) <- rownames(res_5e) <- outcomes 
print(round(res_5d[,-c(5:7)],4))

# Figure with only ITT+controls
# This is the main figure, after reviewer suggestions 
tp <- cbind(res_3d) 
png(file = "Figures/figure-letters-Z2vZ1vZ0-simple-selected.png",
    width = 7, height = 8, units= "in" ,res = 300)
par(mfrow=c(1,1),mar=c(4,10,2,1))
estcol=c(grep("Estimate\\.1",colnames(tp)),grep("Estimate\\.2",colnames(tp)))
secol=c(grep("Std. Error.1",colnames(tp)),grep("Std. Error.2",colnames(tp)))
the.labels=outcome.labels
offsets <- c(.08,-.08)
ys <- -1:-nrow(tp)
xs <- c(-0.1,
        max(.20,max(tp[,estcol]+1.64*tp[,secol])))
plot(xs,c(-.25,-nrow(tp)-0.35)
     ,type="n",bty="n",yaxt="n"
     ,ylab="",xlab="Treatment Effects")	
axis(side=2,at=ys
     ,labels=the.labels 
     ,las=2,tick=F,cex.axis=0.8)
abline(v=0,lty=1,col=gray(0.7))
abline(h=ys,lty=3,col=gray(0.7))
bgs <- c(0,gray(0.5))
pchs <- c(21,24)
cols <- c(1,gray(0.5))
for(i in 1:length(estcol)){
  segments(x0= tp[,estcol[i]]-1.64*tp[,secol[i]],
           x1= tp[,estcol[i]]+1.64*tp[,secol[i]],
           y0=ys+offsets[i],
           y1=ys+offsets[i],
           col= cols[i])
  points(tp[,estcol[i]],ys+offsets[i],bg=bgs[i],pch=pchs[i],cex=0.8,col=1)	
}
legend(x="topright",legend=c("In Person","Instant Messaging"),
       ,pch=c(21,24),col=1,bty="n",pt.bg=c(0,gray(0.5)),lty=1
       ,ncol=1,cex=0.8,inset=-0.025,xpd=NA)
dev.off()

## Summarized plot (with ineligibility and two models only)
## This was the figure in the original submission
## and is now in the supplemental information for the paper
tp <- cbind(res_3d,res_5d) 
png(file = "Figures/figure-letters-Z2vZ1vZ0-ineligibility-selected.png",
    width = 7, height = 8, units= "in" ,res = 300)
par(mfrow=c(1,1),mar=c(4,10,2,1))
estcol=c(grep("Estimate\\.1",colnames(tp)),grep("Estimate\\.2",colnames(tp)))
secol=c(grep("Std. Error.1",colnames(tp)),grep("Std. Error.2",colnames(tp)))
the.labels=outcome.labels
offsets <- c(.13,-.08,.08,-.13)
ys <- -1:-nrow(tp)
xs <- c(-0.1,
        max(.25,max(tp[,estcol]+1.64*tp[,secol])))
plot(xs,c(-.25,-nrow(tp)-0.35)
     ,type="n",bty="n",yaxt="n"
     ,ylab="",xlab="Treatment Effects")	
axis(side=2,at=ys
     ,labels=the.labels 
     ,las=2,tick=F,cex.axis=0.8)
abline(v=0,lty=1,col=gray(0.7))
abline(h=ys,lty=3,col=gray(0.7))
bgs <- c(0,0,gray(0.5),gray(0.5))
pchs <- c(21,24,21,24)
cols <- c(1,1,gray(0.5),gray(0.5))
for(i in 1:length(estcol)){
  segments(x0= tp[,estcol[i]]-1.64*tp[,secol[i]],
           x1= tp[,estcol[i]]+1.64*tp[,secol[i]],
           y0=ys+offsets[i],
           y1=ys+offsets[i],
           col= cols[i])
  points(tp[,estcol[i]],ys+offsets[i],bg=bgs[i],pch=pchs[i],cex=0.8,col=1)	
}
legend(x="top",legend=c("ITT","CACE  ","ITT","CACE  "),title="In Person         Via Messaging"
       ,pch=c(21,24),col=1,bty="n",pt.bg=c(0,0,gray(0.5),gray(0.5)),lty=1
       ,ncol=2,cex=0.8,inset=-0.025,xpd=NA)
dev.off()

## Figure without controlling for ineligibility (in the Supplemental Information)
## Summarized plot (with ineligibility and two models only)
tp <- cbind(res_3a,res_3b,res_5b) 
png(file = "Figures/figure-letters-Z2vZ1vZ0-selected.png",
    width = 7, height = 8, units= "in" ,res = 300)
par(mfrow=c(1,1),mar=c(4,10,2,1))
estcol=c(grep("Estimate\\.1",colnames(tp)),grep("Estimate\\.2",colnames(tp)))
secol=c(grep("Std. Error.1",colnames(tp)),grep("Std. Error.2",colnames(tp)))
the.labels=outcome.labels
offsets <- c(.2,.14,.08,-.08,-.14,-.2)
ys <- -1:-nrow(tp)
xs <- c(-0.1,
        max(.25,max(tp[,estcol]+1.64*tp[,secol])))
plot(xs,c(-.25,-nrow(tp)-0.35)
     ,type="n",bty="n",yaxt="n"
     ,ylab="",xlab="Treatment Effects")	
axis(side=2,at=ys
     ,labels=the.labels 
     ,las=2,tick=F,cex.axis=0.8)
abline(v=0,lty=1,col=gray(0.7))
abline(h=ys,lty=3,col=gray(0.7))
bgs <- c(0,0,0,gray(0.5),gray(0.5),gray(0.5))
pchs <- c(21,24,25,21,24,25)
cols <- c(1,1,1,gray(0.5),gray(0.5),gray(0.5))
for(i in 1:length(estcol)){
  segments(x0= tp[,estcol[i]]-1.64*tp[,secol[i]],
           x1= tp[,estcol[i]]+1.64*tp[,secol[i]],
           y0=ys+offsets[i],
           y1=ys+offsets[i],
           col= cols[i])
  points(tp[,estcol[i]],ys+offsets[i],bg=bgs[i],pch=pchs[i],cex=0.8,col=1)	
}
legend(x="top"
       ,legend=c("ITT","ITT+Ctrls","CACE  ","ITT","ITT+Ctrls","CACE  ")
       ,title="In Person         Via Messaging"
       ,pch=c(21,24,25),col=1,bty="n"
       ,pt.bg=c(0,0,0,gray(0.5),gray(0.5),gray(0.5))
       ,lty=1
       ,ncol=2,cex=0.8,inset=-0.025,xpd=NA)
dev.off()

# Figure with raw difference in means, ITT, and CACE
# This figure is in the supplmental information of the final 
tp <- cbind(res_3a,res_3d,res_5d) 
png(file = "Figures/figure-letters-Z2vZ1vZ0-extended-selected.png",
    width = 7, height = 8, units= "in" ,res = 300)
par(mfrow=c(1,1),mar=c(4,10,2,1))
estcol=c(grep("Estimate\\.1",colnames(tp)),grep("Estimate\\.2",colnames(tp)))
secol=c(grep("Std. Error.1",colnames(tp)),grep("Std. Error.2",colnames(tp)))
the.labels=outcome.labels
offsets <- c(.2,.14,.08,-.08,-.14,-.2)
ys <- -1:-nrow(tp)
xs <- c(-0.1,
        max(.25,max(tp[,estcol]+1.64*tp[,secol])))
plot(xs,c(-.25,-nrow(tp)-0.35)
     ,type="n",bty="n",yaxt="n"
     ,ylab="",xlab="Treatment Effects")	
axis(side=2,at=ys
     ,labels=the.labels 
     ,las=2,tick=F,cex.axis=0.8)
abline(v=0,lty=1,col=gray(0.7))
abline(h=ys,lty=3,col=gray(0.7))
bgs <- c(0,0,0,gray(0.5),gray(0.5),gray(0.5))
pchs <- c(21,24,25,21,24,25)
cols <- c(1,1,1,gray(0.5),gray(0.5),gray(0.5))
for(i in 1:length(estcol)){
  segments(x0= tp[,estcol[i]]-1.64*tp[,secol[i]],
           x1= tp[,estcol[i]]+1.64*tp[,secol[i]],
           y0=ys+offsets[i],
           y1=ys+offsets[i],
           col= cols[i])
  points(tp[,estcol[i]],ys+offsets[i],bg=bgs[i],pch=pchs[i],cex=0.8,col=1)	
}
legend(x="top"
       ,legend=c("ITT","ITT+Ctrls","CACE  ","ITT","ITT+Ctrls","CACE  ")
       ,title="In Person         Via Messaging"
       ,pch=c(21,24,25),col=1,bty="n"
       ,pt.bg=c(0,0,0,gray(0.5),gray(0.5),gray(0.5))
       ,lty=1
       ,ncol=2,cex=0.8,inset=-0.025,xpd=NA)
dev.off()

# Results using the continuous version of the knoweledge index
# In Supplemental information
outcomes <- "know_meifeaturesi_e"
nn <- 1
res_3index <-  bind_rows(lapply(1:nn, FUN = est_itt,fe=F
                                ,ctrls=c(the.ctrls,"ineligible"),the.data=ed)) 
res_5index <- bind_rows(lapply(1:nn, FUN = est_cace,fe=F
                               ,ctrls=c(the.ctrls,"ineligible"),baseline=T,the.data=ed)) 

out <- data.frame(Z1.itt=t(res_3index[,c(1,2,4)]),
                  Z1.cace=t(res_5index[,c(1,2,4)]),
                  Z2.itt=t(res_3index[,c(8,9,11)]),
                  Z2.cace=t(res_5index[,c(8,9,11)]))
xout <- xtable(out)
caption(xout) <- "Treatment Effects on Feature Knowledge Index"
label(xout) <- "tab:index"
align(xout) <- "lrrrr"
print(xout,caption.placement="top")


### Table of coefficients for supplemental materials      ###
### Implementing Benjamini & Hochberg p-value adjustments ###
stack.results <- function(x,model.lab,tt=1,FWERadjust=F){
  to.select <-  which(is.element(names(x),paste(c("Estimate","Std. Error","Pr(>|t|)"),tt,sep=".")))
  x <- x[,to.select]
  if(FWERadjust==T){
    cat("Implementing Benjamini & Hochberg (1995) correction for FWER\n") 
    x[,3] <- p.adjust(x[,3],"BH")
    }
  names(x) <- c("1est","2se","3pval")
  x <- stack(x)
  x$ind <- paste(outcomes,x$ind,sep="")
  x <- x[order(x$ind),]
  names(x)[1] <- model.lab
  return(x)}

the.tab.1 <- merge(
               merge(stack.results(res_3a,"RawITT",tt=1,FWERadjust=T)
                    ,stack.results(res_3d,"ITTwCtrls",tt=1,FWERadjust=T)
                    ,by="ind"),
               stack.results(res_5d,"CACE",tt=1),by="ind")

the.tab.2 <- merge(
              merge(stack.results(res_3a,"RawITT",tt=2,FWERadjust=T)
                  ,stack.results(res_3d,"ITTwCtrls",tt=2,FWERadjust=T)
                  ,by="ind"),
              stack.results(res_5d,"CACE",tt=2),by="ind")

the.tab <- merge(the.tab.1,the.tab.2
                 ,by="ind"
                 ,suffixes=c("InPerson","InstMess"))

xtab <- xtable(the.tab,digits=3)
caption(xtab) <- "Table of Coefficients (with Benjamini & Hochberg adjusted p-values)"
label(xtab) <- "tab:coefs"
print(xtab,include.rownames=F,caption.placement = "top")
